#include "StdAfx.h"
#include "StratumInterp.h"

/*
z = sum(z(i) / (d(i)^k)) / sum(1 / (d(i)^k));
*/

float InterpIDW(const CArray2D<glm::vec3>& cDrillData, UINT uiStratum, float x, float y);

HRESULT StratumInterpIDW(
	const CArray2D<glm::vec3>& cDrillData, 
	CArray3D<glm::vec3>& cStratumData)
{
	glm::vec3 stBoundMin, stBoundMax;
	UINT uiMeshWidth = cStratumData.GetWidth() - 1;
	UINT uiMeshHeight = cStratumData.GetHeight() - 1;

	GetDrillBoundBox(cDrillData, stBoundMin, stBoundMax);

	UINT uiDrillCount = cDrillData.GetHeight();
	UINT uiStratumCount = cDrillData.GetWidth();
	float fCellWidth = (stBoundMax.x - stBoundMin.x) / uiMeshWidth;
	float fCellHeight = (stBoundMax.y - stBoundMin.y) / uiMeshHeight;

	for(UINT uiStratum = 0; uiStratum < uiStratumCount; ++uiStratum)
	{
		for(UINT i = 0; i < uiMeshHeight + 1; ++i)
		{
			for(UINT j = 0; j < uiMeshWidth + 1; ++j)
			{
				cStratumData[uiStratum][i][j].x = stBoundMin.x + j * fCellWidth;
				cStratumData[uiStratum][i][j].y = stBoundMin.y + i * fCellHeight;
				cStratumData[uiStratum][i][j].z = InterpIDW(
					cDrillData, 
					uiStratum, 
					cStratumData[uiStratum][i][j].x, 
					cStratumData[uiStratum][i][j].y);
			}
		}
	}

	return S_OK;
}

float InterpIDW(const CArray2D<glm::vec3>& cDrillData, UINT uiStratum, float x, float y)
{
	UINT uiDrillCount = cDrillData.GetHeight();

	double dNume = 0, dDeno = 0;

	for(UINT uiDrill = 0; uiDrill < uiDrillCount; ++uiDrill)
	{
		float fDist = glm::length(glm::vec2(
			cDrillData[uiDrill][uiStratum].x - x, 
			cDrillData[uiDrill][uiStratum].y - y));
		if(fDist < __EPSILON__)return cDrillData[uiDrill][uiStratum].z;

		dNume += cDrillData[uiDrill][uiStratum].z * 1.0 / (fDist * fDist);
		dDeno += 1.0 / (fDist * fDist);
	}

	return (float)(dNume / dDeno);
}